## ----------------------------------
## Panel A
## ----------------------------------

rm(list = ls())
load("dataBJPOLS_norm.RData")

for(j in c("09", "07", "05")) {
  last.cases <- last.cases.o
  
  date.1 <- paste0("2012-", j, "-06")
  date.2 <- paste0("2016-", j, "-08")
  last.cases$exclude <- (last.cases$start_date_fa_seq > date.1 & last.cases$start_date_fa_seq < "2012-11-06") + (last.cases$start_date_fa_seq > date.2 & last.cases$start_date_fa_seq < "2016-11-08")
  last.cases <- last.cases[over85_flag == 0 & dis_flag == 0 & exclude == 0, ]
  
  inst.1 <- "judgeiv_hd"
  endo.1 <- "pti"
  outc.1 <- "vote_post"
  
  time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)"
  demo.controls <- "age + I(age^2) +  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
  case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"
  
  form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
  form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
  form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))
  
  results <- list()
  run.boot <- T
  
  if(run.boot) {
    set.seed(1231)
    for(i in 1:500) {
      out <- last.cases[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
      last.cases.boot <- last.cases[last.cases$this_id_fa %in% out$V1, ]
      
      m1a1 <- ivreg(form.1, data = last.cases.boot)
      m1a2 <- ivreg(form.2, data = last.cases.boot)
      m1a3 <- ivreg(form.3, data = last.cases.boot)
      
      results[[i]] <- c(m1a1$coefficients["pti"], 
                        m1a2$coefficients["pti"], 
                        m1a3$coefficients["pti"])
    }
    SEs <- round(apply(do.call('rbind', results), 2, sd), 2)
  }
  
  # --------------------------------------------------
  # plot PTI
  # --------------------------------------------------
  # Coefficients
  m1a1 <- ivreg(form.1, data = last.cases)
  m1a2 <- ivreg(form.2, data = last.cases)
  m1a3 <- ivreg(form.3, data = last.cases)
  
  res1_pti <- coefficients(m1a1)[grep("pti", names(coefficients(m1a1)))]
  res2_pti <- coefficients(m1a2)[grep("pti", names(coefficients(m1a2)))]
  res3_pti <- coefficients(m1a3)[grep("pti", names(coefficients(m1a3)))]
  
  m1a1_d <- summary(m1a1, vcov = sandwich, diagnostic = T)
  m1a2_d <- summary(m1a2, vcov = sandwich, diagnostic = T)
  m1a3_d <- summary(m1a3, vcov = sandwich, diagnostic = T)
  
  m1a1_d1 <- m1a1_d$diagnostics[grep("Weak", rownames(m1a1_d$diagnostics)), 3]
  m1a2_d1 <- m1a2_d$diagnostics[grep("Weak", rownames(m1a2_d$diagnostics)), 3]
  m1a3_d1 <- m1a3_d$diagnostics[grep("Weak", rownames(m1a3_d$diagnostics)), 3]
  
  p1 <- data.table(cbind(res1_pti, res2_pti, res3_pti))
  p1$V4 <- 1
  p1$V5 <- 1:nrow(p1)
  p1$V6 <- "Estimate"
  colnames(p1) <- paste0("V", 1:6)
  
  p2 <- data.table(cbind(SEs[1], SEs[2], SEs[3]))
  p2$V4 <- 2
  p2$V5 <- 1:nrow(p2)
  p2$V6 <- "Std. Error"
  
  p3 <- data.table(round(cbind(m1a1_d1, m1a2_d1, m1a3_d1), 2))
  colnames(p3) <- c("V1", "V2", "V3")
  p3$V4 <- 3
  p3$V5 <- 1:nrow(p3)
  p3$V6 <- "F-stat"
  
  table <- rbind(p1, p2, p3)
  table <- table[order(V5, V4), ]
  table$V4 <- table$V5 <- NULL
  p4 <- data.table(m1a1$nobs, m1a2$nobs, m1a3$nobs)    
  p4$V6 <- "N"
  
  table <- rbind(table, p4)
  table <- table[, c("V6", "V1", "V2", "V3")]
  colnames(table) <- c("Variable", "Model 1", "Model 2", "Model 3")
  print(table)
}


## ----------------------------------
## Panel B
## ----------------------------------
rm(list = ls())
load("dataBJPOLS.RData")

last.cases$incap <- NA
last.cases$incap[last.cases$time_period == 1] <- as.numeric(last.cases$pti_incap_12_seq[last.cases$time_period == 1] + last.cases$pci_incap_12[last.cases$time_period == 1] >= 1)
last.cases$incap[last.cases$time_period == 2] <- as.numeric(last.cases$pti_incap_16_seq[last.cases$time_period == 2] + last.cases$pci_incap_16[last.cases$time_period == 2] >= 1)

last.cases <- last.cases[incap == 0, ]

inst.1 <- "judgeiv_hd"
endo.1 <- "pti"
outc.1 <- "vote_post"

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)"
demo.controls <- "age + I(age^2) +  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

results <- list()
boot.run <- T

if(boot.run) {
  set.seed(1231)
  for(i in 1:500) {
    out <- last.cases[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
    last.cases.boot <- last.cases[last.cases$this_id_fa %in% out$V1, ]
    
    m1a1 <- ivreg(form.1, data = last.cases.boot)
    m1a2 <- ivreg(form.2, data = last.cases.boot)
    m1a3 <- ivreg(form.3, data = last.cases.boot)
    
    results[[i]] <- c(m1a1$coefficients["pti"], 
                      m1a2$coefficients["pti"], 
                      m1a3$coefficients["pti"])
  }
  
  SEs <- round(apply(do.call('rbind', results), 2, sd), 2)
}

m1a1 <- ivreg(form.1, data = last.cases)
m1a2 <- ivreg(form.2, data = last.cases)
m1a3 <- ivreg(form.3, data = last.cases)

res1_pti <- coefficients(m1a1)[grep("pti", names(coefficients(m1a1)))]
res2_pti <- coefficients(m1a2)[grep("pti", names(coefficients(m1a2)))]
res3_pti <- coefficients(m1a3)[grep("pti", names(coefficients(m1a3)))]

c1 <- c(res1_pti, SEs[1])
c2 <- c(res2_pti, SEs[2])
c3 <- c(res3_pti, SEs[3])

m1 <- data.frame(rbind(c1,
                       c2,
                       c3))

colnames(m1) <- c("estimate", "std.error")
m1$term <- c("t1", "t2", "t3")
m1_df <- data.frame(as.matrix(m1)) 
m1_df$estimate <- as.numeric(as.character(m1_df$estimate))
m1_df$std.error <- as.numeric(as.character(m1_df$std.error))

m1a1_d <- summary(m1a1, vcov = sandwich, diagnostic = T)
m1a2_d <- summary(m1a2, vcov = sandwich, diagnostic = T)
m1a3_d <- summary(m1a3, vcov = sandwich, diagnostic = T)

m1a1_d1 <- m1a1_d$diagnostics[grep("Weak", rownames(m1a1_d$diagnostics)), 3]
m1a2_d1 <- m1a2_d$diagnostics[grep("Weak", rownames(m1a2_d$diagnostics)), 3]
m1a3_d1 <- m1a3_d$diagnostics[grep("Weak", rownames(m1a3_d$diagnostics)), 3]

p1 <- data.table(cbind(res1_pti, res2_pti, res3_pti))
p1$V4 <- 1
p1$V5 <- 1:nrow(p1)
p1$V6 <- "Estimate"
colnames(p1) <- paste0("V", 1:6)

p2 <- data.table(cbind(SEs[1], SEs[2], SEs[3]))
p2$V4 <- 2
p2$V5 <- 1:nrow(p2)
p2$V6 <- "Std. Error"

p3 <- data.table(round(cbind(m1a1_d1, m1a2_d1, m1a3_d1), 2))
colnames(p3) <- c("V1", "V2", "V3")
p3$V4 <- 3
p3$V5 <- 1:nrow(p3)
p3$V6 <- "F-stat"

table <- rbind(p1, p2, p3)
table <- table[order(V5, V4), ]
table$V4 <- table$V5 <- NULL
p4 <- data.table(m1a1$nobs, m1a2$nobs, m1a3$nobs)    
p4$V6 <- "N"

table <- rbind(table, p4)
table <- table[, c("V6", "V1", "V2", "V3")]
colnames(table) <- c("Variable", "Model 1", "Model 2", "Model 3")
table

## ----------------------------------
## Panel C
## ----------------------------------
rm(list = ls())
load("dataBJPOLS.RData")

last.cases$incap <- NA
last.cases$incap[last.cases$time_period == 1] <- as.numeric(last.cases$pti_incap_12_seq[last.cases$time_period == 1] + last.cases$pci_incap_12[last.cases$time_period == 1] >= 1)
last.cases$incap[last.cases$time_period == 2] <- as.numeric(last.cases$pti_incap_16_seq[last.cases$time_period == 2] + last.cases$pci_incap_16[last.cases$time_period == 2] >= 1)

inst.1 <- "judgeiv_hd * pci_incap_lowp"
endo.1 <- "pti * pci_incap_lowp"
outc.1 <- "vote_post"

time.controls <- "as.factor(court_time1) + as.factor(court_time2) + as.factor(court_dow) + as.factor(severity)"
demo.controls <- "age + I(age^2) +  as.factor(race4) + female + vote_pre + as.factor(noteli) + regis_before"
case.controls <- "as.factor(any_drug) +  as.factor(any_weapon) +  as.factor(any_prop) + as.factor(any_prior_case)"

form.1 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "|", inst.1, "+", time.controls))
form.2 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "|", inst.1, "+", time.controls, "+", demo.controls))
form.3 <- formula(paste(outc.1, "~", endo.1, "+" , time.controls, "+", demo.controls, "+", case.controls, "|", inst.1, "+", time.controls, "+", demo.controls, "+", case.controls))

results <- resultsW <- list()
boot.run <- T

if(boot.run) {
  set.seed(1231)
  for(i in 1:500) {
    out <- last.cases[, this_id_fa[sample.int(.N, .N, TRUE)], by = c("judge_cat", "court_time1")]
    last.cases.boot <- last.cases[last.cases$this_id_fa %in% out$V1, ]
    
    m1a1 <- ivreg(form.1, data = last.cases.boot)
    m1a2 <- ivreg(form.2, data = last.cases.boot)
    m1a3 <- ivreg(form.3, data = last.cases.boot)
    
    results[[i]] <- c(m1a1$coefficients["pti"], 
                      m1a2$coefficients["pti"], 
                      m1a3$coefficients["pti"])
    resultsW[[i]] <- c(m1a1$coefficients["pti:pci_incap_lowp"], 
                       m1a2$coefficients["pti:pci_incap_lowp"], 
                       m1a3$coefficients["pti:pci_incap_lowp"])
  }
  
  SEs <- round(apply(do.call('rbind', results), 2, sd), 2)
  SEsW <- round(apply(do.call('rbind', resultsW), 2, sd), 2)
}

m1a1 <- ivreg(form.1, data = last.cases)
m1a2 <- ivreg(form.2, data = last.cases)
m1a3 <- ivreg(form.3, data = last.cases)

res1_pti <- coefficients(m1a1)[grep("pti", names(coefficients(m1a1)))]
res2_pti <- coefficients(m1a2)[grep("pti", names(coefficients(m1a2)))]
res3_pti <- coefficients(m1a3)[grep("pti", names(coefficients(m1a3)))]

c1 <- c(res1_pti, SEs[1])
c2 <- c(res2_pti, SEs[2])
c3 <- c(res3_pti, SEs[3])

m1 <- data.frame(rbind(c1,
                       c2,
                       c3))

colnames(m1) <- c("estimate", "std.error")
m1$term <- c("t1", "t2", "t3")
m1_df <- data.frame(as.matrix(m1)) 
m1_df$estimate <- as.numeric(as.character(m1_df$estimate))
m1_df$std.error <- as.numeric(as.character(m1_df$std.error))

m1a1_d <- summary(m1a1, vcov = sandwich, diagnostic = T)
m1a2_d <- summary(m1a2, vcov = sandwich, diagnostic = T)
m1a3_d <- summary(m1a3, vcov = sandwich, diagnostic = T)

m1a1_d1 <- m1a1_d$diagnostics[grep("Weak", rownames(m1a1_d$diagnostics)), 3]
m1a2_d1 <- m1a2_d$diagnostics[grep("Weak", rownames(m1a2_d$diagnostics)), 3]
m1a3_d1 <- m1a3_d$diagnostics[grep("Weak", rownames(m1a3_d$diagnostics)), 3]

p1 <- data.table(cbind(res1_pti, res2_pti, res3_pti))
p1$V4 <- 1
p1$V5 <- 1:nrow(p1)
p1$V6 <- "Estimate"
colnames(p1) <- paste0("V", 1:6)

p2 <- data.table(rbind(SEs, SEsW))
colnames(p2) <- c("V1", "V2", "V3")
p2$V4 <- 2
p2$V5 <- 1:nrow(p2)
p2$V6 <- "Std. Error"

p3 <- data.table(round(cbind(m1a1_d1, m1a2_d1, m1a3_d1), 2))
colnames(p3) <- c("V1", "V2", "V3")
p3$V4 <- 3
p3$V5 <- 1:nrow(p3)
p3$V6 <- "F-stat"

table <- rbind(p1, p2, p3)
table <- table[order(V5, V4), ]
table$V4 <- table$V5 <- NULL
p4 <- data.table(m1a1$nobs, m1a2$nobs, m1a3$nobs)    
p4$V6 <- "N"

table <- rbind(table, p4)
table <- table[, c("V6", "V1", "V2", "V3")]
colnames(table) <- c("Variable", "Model 1", "Model 2", "Model 3")
table